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Abstract. We report a current status of our radiation-magnetohydrodynamic code for the 
study of core-collapse supernovae. In this contribution, we discuss the accuracy of our newly 
developed numerical code by presenting the test problem in a static background model. We 
also present the application to the spherically symmetric core-collapse simulations. Since close 
comparison with the previously published models is made, we are now applying it for the study 
of magnetorotational core-collapse supernovae. 



1. Introduction 

Recently, multidimensional studies and simulations of core-collapse supernovae have come into 
blossom again since 1990's, when the direct observations of global asymmetry in SN 1987A were 
reported. This trend may be partly because spherically symmetric supernova simulations have 
not yet produced explosions (e.g., pQ) albeit with the sophistication of the neutrino-transport 
method including the state-of-the-art nuclear physics, and the detailed weak interaction rates. 
Meanwhile, there have been the accumulating observations after SN 1987A implying that the 
core-collapse supernovae are generally aspherical (e.g., pj). The degree of asymmetry also rises 
as a function of time, which has been interpreted as an evidence that the inner portion of the 
explosion are strongly aspherical. 

A leap beyond the spherical models seems indeed meaningful, because asphericities in the 
supernova core should have influence on the nucleosynthesis, neutrino and gravitational-wave 
emissions 3 . Before the advent of the observations using neutrinos and gravitational waves as 
new eyes, one hopes to clarify the explosion mechanism of aspherical supernovae. 

So far, many physical ingredients to produce aspherical explosion have been investigated, 
such as convection, possible density inhomogeneities formed prior to core-collapse, rotation and 
magnetic fields |3] (see jH] for a collective reference), more recently standing shock instability 
[5] and the excitations of g-modes in the protoneutron star Whatsoever the origin of the 
asphericity, the neutrino heating mechanism should play a key role to drive explosions. Thus 
multidimensional neutrino transport simulations are indispensable. In this contribution, we 
report the current status of our radiation transport calculations. 



2. Basic equations of radiation-magnetohydrodynamics 

Magnetohydrodynamic (MHD) equations including the neutrino-matter coupling read, 
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The left-hand side of the above equations is MHD evolution equations, which we numerically 
integrate by the Newtonian explicit code ZEUS2D 7 . Equation (^Q) is the Euler equation, where 
p is the mass density, d/dt is the Lagrange derivative, v is the fluid velocity, p is the matter 
pressure, B is the magnetic fields solved by the induction equation : dB/dt = V x (v x B) 
keeping the field divergence- free (V-B = 0) using the CT scheme, $ is the gravitational potential 
including the general relativistic correction. Equation @ and (j3J) expresses the evolution of 
specific internal energy (e) and electron fraction (Y e ), respectively. As for the equation of state, 
we incorporate a tabulated one based on relativistic mean field theory |S]. 

The right-hand side of the above equations expresses the transfer of momentum (Sm), 
exchange of energy (Se) from neutrinos to matter, and the change of electric charge (Sy e )- As for 
the neutrino reactions, we include the so-called standard set denoted in table 1 in plus nucleon 
bremsstrahlung ^Uj. Corrections to the standard neutrino opacities, such as the detailed reaction 
kinematics of nucleon thermal motions, recoil, and weak magnetism, pair-annihilation/creations 
among different neutrino flavors, the quenching of the coupling constant above nuclear density, 
and the modified electron capture rates to the nuclei in the excitation states and so on (see (T^ ) , 
we omit them for simplicity in the simulation presented here. To determine the right-hand side 
of the above equations, one has to know the neutrino distribution function f(t,r,p) by solving 
the Boltzmann equation with the collisional terms including the above neutrino reactions in its 
right-hand side. Fully angle-dependent solutions in multidimensional models (even in 2D) are 
computationally prohibitive, we employ the flux limited diffusion approximation to relate the 
first and zeroth angular moment of the neutrino distribution function, 

V>i = ( ^ ) = -A • [Wo " Mik ~ Ci] , (4) 

where ipQ and xj)\ is the zeroth and first angular moment of the neutrino distribution function, 
respectively. See equation (A25) in jUj for the definition of A\ and C\. Here A is a flux limiter 
generalized to 2D as follows, 

A r 
Ao 
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whose components are given by A r = (3A*)/(3 + A*|V r '0o|/V ; o)? and Ag = (3A*)/(3 + A* | V^V'ol/V'o) 
where A* is the transport mean free path (see Equation (A26) in V r and V# represents d/dr 
and d/(rd0), respectively. With equation (JIJ), one can determine Vo by solving the following 
zeroth angular moment equation of Boltzmann equation, 

+ V • Vi + -V • v u ^-Vo = i(l - ifo) - t£ + A^o + B ■ if,! + C , (6) 
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which is often referred to as a multi-group flux limited diffusion (MGFLD) equation for neutrinos 
(again see appendix A in (Hj for the definition of unmentioned variables). For the time evolution 



of transport, we solve the neutrino-matter coupling (equations (J2J (the pV • v term is treated 
in an operator-splitting fashion), ([5]). and © implicitly with performing the Newton Raphson 
iteration until 5ipo, 5Y e , and 5(e/p) converges to a certain value. More details for the numerical 
implementation of our code will be presented in the forthcoming paper. 

3. Static background 

Before we apply the newly developed code to the magnetized and rotating supernova simulations, 
we need to verify the accuracy by the test calculations. Taking the profiles of the density, 
temperature, and the electron fraction from Bruenn's ID model after core bounce |12| as a 
background, we calculate the neutrino distribution functions of each species (f e , P e , u x , v x ) by 
our 2D computations and compare them with the ones obtained in ID simulation [12] ■ We use 
the 16 energy mesh points which is logarithmically uniform and covers 0.95 - 255 MeV. In Figure 
n the distribution functions of v e and v x at each energy bin as a function of radius is shown. 
One can see the good agreement with each other, except for the small discrepancies at ~ 20 km 
for the lower energy bins. This is also the case for the electron neutrinos, which is not presented 
here. This may be mainly due to the difference of the numerical implementation of the neutrino 
transport, and of the employed equation of state. Since close comparison is made in the static 
background computation, we move on to apply the code to dynamical models. 
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Figure 1. Anti-electron (left) and mu-tau (right) neutrino distribution functions at each energy 
bin as a function of the stellar radius (from to 80 km, equidistantly plotted). Note in the figure 
that = 1, 16) represents the neutrino energy group and that "ID" represents the results by 



4. Spherically symmetric simulation 

Here we examine the accuracy of our code by the spherically symmetric core-collapse simulations 
of 15 Mq progenitor star (Woosley & Weaver 1995). In Figured the velocity (top right), entropy 
(bottom left), Y e (bottom right) distributions at 10 msec after bounce with the evolution of Y e up 
to ~ 40 msec after bounce (top left), are shown. The global profiles are well in good agreement 
with the previous spherically symmetric calculations pQ. Discrepancies with other calculations 
could be originated from the differences of employed microphysics, numerical treatment of 
neutrino transport, general relativistic effects, which seem rather common from models to models 
PQ. More details about the tests of our newly developed code will be presented soon in elsewhere, 
with its applications to the magnetorotational core-collapse simulations of massive stars. 




Figure 2. Top left panel is the time evolution of Y e up to ~40 msec after bounce (small blank 
is due to the data loss). Top right, bottom left, bottom right shows the velocity, entropy, Y e 
distribution at 10 msec after bounce. Compare panels (a), (c), (d) in Figure 20 in [T3] . 



5. Outlook 

Although 0(v/c) terms in the transport equations, which are one of the most important 
ingredients for determining the failure or success of the explosion are neglected in the 
current simulations, we are now implementing them, simultaneously running the MHD collapse 
models without to see the difference. It is better to employ more detailed neutrino microphysics. 
Steady progresses with the parallelization of the code are made to test the predictions in our 
earlier works. 
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